rm(list=ls())
library(foreign)
library(xtable)
library(doBy)




###REPLACE WITH YOUR FILE PATH HERE using setwd("path")

##################################
##################################
##################################
##### Berrebi and Klor (2008)
##################################
##################################
##################################


d<-read.dta("bk.dta")

###Replicate Table 5 Model 1 in Berrebi and Klor (2008)
##Authors' Stata code: xtreg  area_rel_right_bloc tot_fat_3months, fe cluster(area_code)

reg<-lm(area_rel_right_bloc ~tot_fat_3months  + factor(area_code) , data= d)
summary(reg)
coef(reg)[1:5]##replicates the 0.45 pp effect
dim(reg$model)

x<-"tot_fat_3months"

##store coefficient on treatment
b<-coef(reg)[x]
b
##get estimation data
est<-reg$model
dim(est)

table(est[,x])

##get SD of X 
sdx<-sd(est[,x])

##residualize with respect to fixed effects
x.resid<-lm(tot_fat_3months ~ est[,"factor(area_code)"], data=est)$residuals
sdxresid<-sd(x.resid)

#the counterfactual discussed in the article (a one-unit shift)
count.obs<-1 ##one-unit shift

#ratio of discussed counterfactual to our preferred counterfactual
rat.count<-count.obs/(sdxresid)
rat.count


results<-as.character(c("Berrebi and Klor (2008)",round(b,3), round(sdx,3), round(sdxresid,3),   round(((sdxresid-sdx)/(sdx))*100,3 ), round(rat.count,3)   ))


results1<-as.data.frame(t(results))
names(results1)<-c("study","beta", "sdx","sdxresid","%difference","counter1/counter2")
res2<-results1
save(res2, file="comparisontable.Rdata")

res2







##################################
##################################
##################################
#####Gasper and Reeves (2011)
##################################
##################################
##################################


d <- read.csv("gr.csv")

head(d)

##Replicate Table 1, Column 3 in Gasper and Reeves (2011)\\\
reg <- lm(inc.twpct ~ presgov.twpct + govvote.lag +  propertypercapl.6mo + disdecs.all.6mo + turndowns.6mo + medincK + as.factor(countyfips) + as.factor(year), data = d)
summary(reg)$coefficients[1:10,]

table(d$turndowns.6mo)

x<-"turndowns.6mo"
b<-coef(reg)[x]
b

##get estimation data
est<-reg$model
table(est[,x])


##residualize with respect to fixed effects
x.resid<-lm(turndowns.6mo~est[,"as.factor(countyfips)"] + est[,"as.factor(year)"], data=est)$residuals
sdxresid<-sd(x.resid)
sdx<-sd(est$turndowns.6mo)

#the counterfactual discussed in the article (one unit shift)
count.obs<-1 #one-unit shift

#ratio of discussed counterfactual to our preferred counterfactual
rat.count<-count.obs/(sdxresid)
rat.count

results<-as.character(c("Gasper and Reeves (2011)",round(b,3), round(sdx,3), round(sdxresid,3),   round(((sdxresid-sdx)/(sdx))*100,3 ), round(rat.count,3)  ))

results1<-as.data.frame(t(results))
names(results1)<-c("study","beta", "sdx","sdxresid","%difference","counter1/counter2")

results1

load("comparisontable.Rdata")

 res2<-rbind(res2, results1)
 save(res2, file="comparisontable.Rdata")
res2


##################################
##################################
##################################
#####Ichino and Nathan (2013)
##################################
##################################
##################################



d <- read.csv("in.csv", header=TRUE)
dim(d)
head(d)

d <- d[,c("npp2008ps_pres_p",  "ndc2008ps_pres_p", "akan_p_poly", "ewe_p_poly", "moledagbon_p_poly", "public_semipublic_p", "dev_factor2", "akan_30km_l_p", "akan_20km_l_p", "akan_40km_l_p", "ewe_30km_l_p", "ewe_20km_l_p", "ewe_40km_l_p", "h_30rad_e", "c230_id_h", "otherethn_p_poly", "totalvotes2008ps_pres", "ethfrac_30km_l", "akan_p")]
dim(d)
d <- na.omit(d)
dim(d)


#Replicate Table 2, Model 2 of Ichino and Nathan ####
reg <- lm(npp2008ps_pres_p ~ akan_p_poly + moledagbon_p_poly + otherethn_p_poly + public_semipublic_p + dev_factor2 + akan_30km_l_p + factor(c230_id_h), data=d, weights= totalvotes2008ps_pres)

est<-reg$model
x<-"akan_30km_l_p"
b<-coef(reg)[x]
b


#use weights since the used them in the model (though it makes very little difference)
x.resid<-lm(akan_30km_l_p ~factor(c230_id_h), data=d, weights= totalvotes2008ps_pres)$residuals


sdx<-sd(est[,x])

sdxresid<-sd(x.resid)
#the counterfactual discussed in the article
count.obs<-sdx #one-SD shift
count.obs


#ratio of discussed counterfactual to our preferred counterfactual
rat.count<-count.obs/(sdxresid)
rat.count


results<-as.character(c("Ichino and Nathan (2013)",round(b,3), round(sdx,3), round(sdxresid,3),   round(((sdxresid-sdx)/(sdx))*100,3 ), round(rat.count,3)  ))


results1<-as.data.frame(t(results))
names(results1)<-c("study","beta", "sdx","sdxresid","%difference","counter1/counter2")
results1

load("comparisontable.Rdata")

 res2<-rbind(res2, results1)
 save(res2, file="comparisontable.Rdata")
res2




##################################
##################################
##################################
##### Scheve and Stasavage (2012)
##################################
##################################
##################################


d <- read.dta(file="ss.dta")

##Authors' Stata code for model: xi: reg  firsttopitaxrate2 himobpopyear2pl1 unisuffragel1 leftexec2l1 rgdppcl1 i.ccode i.hdecadec,cluster(ccode)


#Replicate Table 2 model 1 in Scheve and Stasavage (2012)
reg <- lm(firsttopitaxrate2 ~ himobpopyear2pl1 + unisuffragel1 + factor(ccode) + factor(hdecadec) , data = d)

##save estimation data
est<-reg$model


#treatment
x<-"himobpopyear2pl1"

#how many levels to the treatment?
table(est[,x], exclude=NULL)


dim(est)
head(est)

b <- reg$coefficients[x]
b

table(est[,x])

x.resid <- lm(himobpopyear2pl1 ~ est[,"factor(ccode)"] + est[,"factor(hdecadec)"] , data = est)$residuals

sdx<-sd(est$himobpopyear2pl1)
sdxresid<-sd(x.resid)


##what share of observations experience the max shift
table(est[,x])["1"]/nrow(est) ##about 1.5 %

##shift of 1 represents a 5 year war mobilization within country, after applying year dummies. thats about the 97th percent of the disribution of within country ranges

#the counterfactual discussed in the article
count.obs<-1 #(one-unit/min to max shift)

#ratio of discussed counterfactual to our preferred counterfactual

rat.count<-count.obs/(sdxresid)
rat.count


results<-as.character(c("Scheve and Stasavage (2012)",round(b,3), round(sdx,3), round(sdxresid,3),   round(((sdxresid-sdx)/(sdx))*100,3 ), round(rat.count,3)  ))


results1<-as.data.frame(t(results))
names(results1)<-c("study","beta", "sdx","sdxresid","%difference","counter1/counter2")
results1

results1

load("comparisontable.Rdata")

 res2<-rbind(res2, results1)
 
 save(res2, file="comparisontable.Rdata")
res2



##################################
##################################
##################################
#####Snyder & Stromberg (2010)
##################################
##################################
##################################



load("ss2.Rdata")
d<-congrecall
dim(d)

##Replicate Table 4, Row 2, Column 1 in Snyder & stromberg (2010)
reg<-lm(cong_recall ~ cov2c + factor(year) + factor(incnum), data=d)

est<-reg$model
x<-"cov2c"
b<-coef(reg)[x]
b
sdx<-sd(est$cov2c)
(range(est$cov2c)[2]-range(est$cov2c)[1])*b

##Calculate Mean-deviated X with respect to years

##Residualize with respect to fixed effects
x.resid<-(lm(cov2c~est[,"factor(year)"] + est[,"factor(incnum)"], data=est))$residuals
sdxresid<-sd(x.resid)


est$incnum<-est[,"factor(incnum)"]
est$year<-est[,"factor(year)"]
head(est)

est2<-summaryBy(cov2c ~incnum, data=est,FUN=range)
est2$range<-est2[,3] - est2[,2]
summary(est2$range)

##shift the size of mean of within-incumbent ranges
round(mean(est2$range),3 )
b*round(mean(est2$range),3 )


##shift the size of 95th percentile of within-incumbent ranges
quantile(est2$range, probs=.95)
round(b*quantile(est2$range, probs=.95),3 )


#the counterfactual discussed in the article
count.obs<-1 #one-unit/min to max shift

#ratio of discussed counterfactual to our preferred counterfactual
rat.count<-count.obs/(sdxresid)
rat.count


results<-as.character(c("Snyder and Stromberg (2010)",round(b,3), round(sdx,3), round(sdxresid,3),   round(((sdxresid-sdx)/(sdx))*100,3 ), round(rat.count,3)   ))


results1<-as.data.frame(t(results))
names(results1)<-c("study","beta", "sdx","sdxresid","%difference","counter1/counter2")
results1


load("comparisontable.Rdata")

res2<-rbind(res2, results1)
 



save(res2, file="comparisontable.Rdata")
res2


##store necessary columns
res2<-res2[,c("study","beta","sdx","sdxresid","%difference","counter1/counter2")]
res2
 save(res2, file="comparisontable.Rdata")
xtable(res2)




##generate Snyder and Stromberg Figure for paper (Figure 1)

res<-as.data.frame(cbind.data.frame(v1=na.omit(est$cov2c-mean(est$cov2c, na.rm=T)),v1resid=x.resid-mean(x.resid, na.rm=T)  ))
head(res)

#Variation in the Treatment in Snyder and Stromberg (2010)\n Before and After Fixed Effect
##centering both distributions at zero
pdf(file="snyder_strom_density.pdf", height=6, width=8)
par(mfrow=c(1,2))
plot(density(x= res$v1, bw=.01),  col="black", cex=.3, lty=2, main="Treatment Distribution", xlab="Treatment (Centered on Zero)", ylim=c(0,25), axes=F)
lines(density(x= res$v1resid, bw=.01), col="red", cex=.7, xlab="")
text(x=.22, y=22, labels="Distribution\nafter Fixed\nEffects", col="red")
text(x=.22,y=19.5, labels=paste("SD=",round(sd(res$v1resid),3),sep=""), cex=1, col="red" )
text(x=.3, y=5.5, labels="Original\nDistribution", col="black")
text(x=.3,y=3.5, labels=paste("SD=",round(sd(res$v1),3),sep=""), cex=1 )
axis(1, seq(-1,1, by=.2))
axis(2, seq(0,40, by=10), las=2)

hist(est2$range,col="grey", main="Within-Unit Ranges\nof Treatment",xlab="Within-Incumbent Ranges", ylim=c(0,1000), xlim=c(0,1),axes=F)
axis(1, seq(0,1, by=.2))
axis(2, seq(0,1000, by=200), las=2)
abline(v=mean(est2$range), lty=2)
abline(v=quantile(est2$range, probs=.95), lty=2)
abline(v=1, lty=2)
text(x=.2,y=900,labels="mean")
text(x=.7,y=900,labels="95th\npercentile")
text(x=.68,y=500,labels="counterfactual\ndiscussed in text")
arrows(x0=.18,x1=mean(est2$range), y0=875, y1=875, length=.1)
arrows(x0=.7,x1=quantile(est2$range, probs=.95), y0=840, y1=840, length=.1)
arrows(x0=.7,x1=1, y0=440, y1=440, length=.1)
dev.off()


##what share of ranges are zero?
mean(I(est2$range==0))##44%




